#include "dna.h"
#include "getargs.h"
#include "libheader.h"

struct IndexStru{
  int dna_no;
  short position;
  unsigned long code;
  IndexStru* next;
};

unsigned long encode1(const unsigned short* seq, int k_length){
  unsigned long value1=0x0;
  
  for(int i=0;i<k_length;i++){
    value1 = value1 << 2;
    value1 = value1 | (unsigned long)seq[i];
  }
  
  return value1;
}

unsigned long encode2(const unsigned short* seq, int k_length){
  unsigned long value2=0x0;
  
  for(int i=k_length-1;i>=0;i--){
    value2 = value2 << 2;
    value2 = value2 | (unsigned long)(0x3 - seq[i]);
  }
  
  return value2;
}

int getRead(vector<DNA*> &read_vec, const char* in_filepath){
  ifstream read_file(in_filepath);
  if (!read_file) {
    cout << "prot file not found : " << in_filepath<<endl;
    return 1;
  }
   
  const int buffer_len = 67108864;
  char* line;
  char* line_buf=new char[buffer_len];
  char* seq_buf=new char[buffer_len];
  char* seq_name=new char[buffer_len];
  seq_name[0]=0;
  string tmp_seq("");
  
  DNA *dna = NULL;
  while (!read_file.eof()) {
    read_file.getline(line_buf,buffer_len);
    
    int start=0;
    while(isspace(line_buf[start]) && start < (int)strlen(line_buf)){
      start++;
    }
    
    if(start==(int)strlen(line_buf))
      continue;
    else
      line=&line_buf[start];
    
    if (line[0] == '>') {
      //read a new sequence
      int seq_len = tmp_seq.length();
      if ( seq_len > 0 ) {
	dna = new DNA( strdup(seq_name), strdup(tmp_seq.c_str()) );
	read_vec.push_back(dna);
	tmp_seq = "";
      }
      strncpy(seq_name, &line[1], strlen(line)-1);
    } else {
      int end = buffer_len - start;
      int current = 0;
      
      while (current < end && line[current] != 0) {
	if ( line[current] == '#' )
	  break; // ignore comment
	char c = toupper(line[current]);
	if ( 'A' <= c && c <= 'Z' ) {
	  seq_buf[current]=c;
	  current ++;
	} else {
	  cout<<"In main.cpp: incorrect character \""<<line[current]<<"\""<<endl;
	  break;
	}
	
      }
      seq_buf[current]=0;
      tmp_seq += seq_buf;
    }
  }

  // last sequence
  dna = new DNA( strdup(seq_name), strdup(tmp_seq.c_str()) );
  read_vec.push_back(dna);
  tmp_seq = "";
  //quence seq(seq_name, tmp_seq.c_str() );
  //quences.push_back(seq);
  //q_file.close();

  read_file.close();
  
  delete[] line_buf;
  delete[] seq_buf;
  delete[] seq_name;
  
  return 0;
}

//get file name in a directory
int getdir (string dir, vector<string> &files)
{
    DIR *dp;
    struct dirent *dirp;
    if((dp  = opendir(dir.c_str())) == NULL) {
        cout << "Error(" << errno << ") opening " << dir << endl;
        return errno;
    }

    while ((dirp = readdir(dp)) != NULL) {
        files.push_back(string(dirp->d_name));
    }
    closedir(dp);
    return 0;
}

//search motif in protein complex with 3D structure
int main(int argc,const char** argv){
  time_t start_time = time(0);
  
  srand(time(0));
  
  struct ParaFormat para_format[] =
  {
    { 0, "in_filepath", "input file", "t", 1, 's'},
    { 1, "out_filepath", "output file name (txt)", "O", 1, 's'},
  };
  
  Getargs args_list(argc,argv,para_format,2);
  
  if(args_list.isError()){
    args_list.print_list(cout);
    return -1;
  }
  
  const char* out_filepath = NULL;
  const char* in_filepath = NULL;
  
/*****************************************************************************/
/*                                                                           */
/* Read in parameters                                                        */
/*                                                                           */
/*****************************************************************************/
  for(int i=0;i<(int)args_list.getParaValSize();i++){ 
    switch((args_list.getPara(i))->getIndex()){
      
      //case0: input file path
      case 0:
        in_filepath=(args_list.getPara(i))->getStrVal();
        break;
      
      //case 1: output file path
      case 1:
        out_filepath=(args_list.getPara(i))->getStrVal();
        break;
    }    
  }
  
  if(in_filepath==NULL){
    cout<<endl;
    cout<<"Error! Please input the path for the read sequences (fasta) by the -t command."<<endl;
    cout<<"You can use \""<<argv[0]<<" -h\" to get the list of possible input command."<<endl;
    cout<<endl;
    return -1;
  }

  if(out_filepath==NULL){
    cout<<endl;
    cout<<"Error! Please input the path for the output file (txt) by the -O command."<<endl;
    cout<<"You can use \""<<argv[0]<<" -h\" to get the list of possible input command."<<endl;
    cout<<endl;
    return -1;
  }
  
/*****************************************************************************/
/*                                                                           */
/* Read in read sequences                                                    */
/*                                                                           */
/*****************************************************************************/

  vector<DNA*> in_read_vec;
  getRead(in_read_vec, in_filepath);

/*  
  for(int i=0;i<(int)read_vec.size();i++){
    cout<<'>'<<read_vec.at(i)->getName()<<endl;
    cout<<read_vec.at(i)->getSeq()<<endl;
  }
*/
  int k_length=32;
  cout<<"Building hash table for "<<in_filepath<<endl; 
  cout<<"k-mer length = "<<k_length<<endl;
  
  unsigned int hash_table_size = 268435456; //0.25 M
  IndexStru** hash_table=new IndexStru*[hash_table_size];
  for(int i=0;i<(int)hash_table_size;i++){
    hash_table[i]=NULL;
  }
/*  
  unsigned short* tmp=new unsigned short[4];
  tmp[0]=tmp[1]=0;
  tmp[2]=tmp[3]=1;
  cout<<encode(tmp,4)<<endl;
*/  
  for(int i=0;i<(int)in_read_vec.size();i++){
    const unsigned short* seq_short=in_read_vec.at(i)->getShortRep();
    int end=in_read_vec.at(i)->getSeqLength()-k_length;
    for(int j=0;j<=end;j++){
      unsigned long code1=encode1(&seq_short[j], k_length);
      unsigned long code2=encode2(&seq_short[j], k_length);
      unsigned long code = (code1>code2?code1:code2);
      unsigned long hash_val=(code1^code2)%hash_table_size;
      IndexStru* new_entry=new IndexStru();
      new_entry->dna_no=i;
      new_entry->position=j;
      new_entry->code=code;
      new_entry->next=NULL;
      
      if(hash_table[hash_val]==NULL){
        hash_table[hash_val]=new_entry;
      } else {
        IndexStru* cur=hash_table[hash_val];
        while(cur->next!=NULL){
          cur=cur->next;
        }
        cur->next=new_entry;
      }
    }    
  }
  
  cout<<"finished building hash table"<<endl; 
  
  //for each other file, compare the k-mers
  int pos=strlen(in_filepath)-1;
  //cout<<in_filepath<<endl;
  while(in_filepath[pos]!='/' && pos>0){
    pos--;
  }
  //cout<<"pos = "<<pos<<endl;
  string dir = string(in_filepath, pos);
  //cout<<dir<<endl;
  
  vector<string> files_vec = vector<string>();
  getdir(dir,files_vec);

  for (unsigned int file_no = 0;file_no < files_vec.size();file_no++) {
    //cout << files_vec[file_no] << endl;
    
    //cout<<files_vec[file_no]<<endl;
    if (files_vec[file_no].length()>3 && files_vec[file_no].compare(files_vec[file_no].length()-3,3,".fa") == 0){
      string file_path=string(dir.c_str());
      file_path+='/';
      file_path+=files_vec[file_no];
      
      cout<<"processing "<<file_path<<endl;
      
      vector<DNA*> comp_read_vec;
      getRead(comp_read_vec, file_path.c_str());
      
      for(int i=0;i<(int)comp_read_vec.size();i++){
        const unsigned short* seq_short=comp_read_vec.at(i)->getShortRep();
        int end=comp_read_vec.at(i)->getSeqLength()-k_length;
        for(int j=0;j<=end;j++){
          unsigned long code1=encode1(&seq_short[j], k_length);
          unsigned long code2=encode2(&seq_short[j], k_length);
          unsigned long code = (code1>code2?code1:code2);

          unsigned long hash_val=(code1^code2)%hash_table_size;
          
          IndexStru* cur=hash_table[hash_val];
          while(cur!=NULL){
            if(cur->code == code){
              bool is_print=true;
              for(int k=0;k<k_length;k++){
                if(in_read_vec[cur->dna_no]->getSeq()[cur->position+k]=='N'
                   || comp_read_vec.at(i)->getSeq()[j+k] == 'N'){
                  is_print=false;
                  break;   
                }
              }
              
              if(is_print){
                cout<<in_read_vec[cur->dna_no]->getName()<<endl;
                
                cout<<"DNA "<<cur->dna_no<<" position "<<cur->position <<"\tmatchs"<<endl;
                cout<<comp_read_vec[i]->getName()<<endl;
                cout<<files_vec[file_no]<<"\tDNA "<<i<<" position "<<j<<endl;
                
                cout<<cur->code<<"\t";
                for(int k=0;k<k_length;k++){
                  cout<<in_read_vec[cur->dna_no]->getSeq()[cur->position+k];
                }
                cout<<endl;
                cout<<code<<"\t";
                for(int k=0;k<k_length;k++){
                  cout<<comp_read_vec.at(i)->getSeq()[j+k];
                }
                cout<<endl;
                cout<<endl;
              }
            }
            cur=cur->next;
          }
        }    
      }
    }
    
  }  
  
  cout<<"total running time is : "<<(time(0) - start_time)/3600 <<" hr"<<endl;
  
  return 0;

}
